A new Godunov-type numerical code with a 
non-linear Riemann solver for equations of 
relativistic hydrodynamics 

Pavlo V. Tytarenko, 

Astronomical Observatory of Kyiv Shevchenko University, 3 Observatorna St., 

Kiev, 04053, Ukraine 

Iurii A. Karpenko, Yury M. Sinyukov 

Bogolyubov Institute for Theoretical Physics, lJ^-b Metrolohichna St., Kiev, 03680, 

Ukraine 



Abstract 

We present a second-order upwind numerical scheme for equations of relativistic 
hydrodynamics with a source term. A new non-linear Riemann solver is constructed. 
Solution of a Riemann problem on a cells boundary is based on exact relations in 
case of zero tangential velocities. In this sense our solver is "exact". In case of non- 
zero tangential velocities a reasonable approximation is made that allows to avoid 
solution of very complicated exact equations. The scheme is tested on several one- 
and three-dimensional solutions demonstrating a good performance for both strong 
shocks and strong rarefactions. 
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1 Introduction 

Hydrodynamic flows with relativistic velocities are essential feature of modern 
astrophysics and physics of ultra-relativistic mucleus-nucleus collisions. De- 
scription of such objects as active galactic nuclei, quasars and micro quasars, 
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supernova collapse and explosion, gamma ray bursts, pulsar winds, extragalac- 
tic jets is impossible without using of relativistic hydrodynamics (see e.g. Refs. 
[2l)j45f28j ). Also relativistic hydrodynamics is widely used for description of 
multiparticle production in high energy nucleus-nucleon collisions in the way 
firstly proposed by Landau [27]. Since then 3D hydrodynamics is applied to 
study the hadronic spectra and momentum correlations, radial and elliptic 
flows in relativistic processes of nuclear collisions (see, e.g. review by Kolb & 
Heinz [25]). Diversity of complex micro- and astro- physical relativistic pro- 
cesses requires development of reliable and universal (i.e. applicable to all 
range of initial thermodynamic conditions and velocity distributions) numer- 
ical methods for equations of relativistic hydrodynamics 

Relativistic hydrodynamics numerical calculations have already a rich history. 
The first methods of numerical treatment of equations of relativistic hydrody- 
namics were based on introduction into equations an artificial viscosity |44j . 
This allows to make initially hyperbolic equations of hydrodynamics parabolic 
and as a result discontinuous features like shock waves become continuous. But 
in such methods viscosity should be chosen separately for every particular task 
introducing thus an element of art in science and making corresponding code 
difficult for application. Another difficulties are loosing energy due to viscosity 
and problems with ultrarelativistic regime (see Ref. e.g. [1]). And finally there 
is a danger of losing causality when one passes to the relativistic equations of 
parabolic type [361124"] . 

Eventually artificial viscosity methods were replaced by more advanced ones. 
One of well-known and successful methods is SPH (Smooth Particle Hydro- 
dynamics), see Refs. [30|16] . The obvious disadvantages of SPH are the use of 
Lagrangian coordinates, problems with resolution at low densities and and still 
need for an artificial viscosity at high Mach numbers. Another class is HRSC 
(High Resolution Shock Capturing) schemes based on Godunov method. Our 
method is an HRSC one. The main idea of Godunov method is introduction of 
Riemann problems - ones of the simplest possible initial value problems for hy- 
perbolic systems - at cell boundaries in order to compute fluxes. But complete 
solution of the Riemann problem is a non trivial task. Thus many methods 
have been developed in order to avoid solution of the Riemann problem. A 
comprehensive review can be found in Ref. [32] . 

The most obvious way of simplifying the task is a proper linearization of ini- 
tially non-linear system. Roe-type solvers (following the original work by Roe 
[39] ) reduce the nonlinear system of hydrodynamics equations to a locally lin- 
ear system where the linear Jacobian matrix is chosen in such a way that a 
non-linear shock is a single discontinuity for the linear system too. General- 
ization of the method to equations of relativistic hydrodynamics was made in 
Ref. [13J. Another method of local linearization based on using primitive vari- 
ables instead of conservative ones and applying either HLLE-type or Roe-type 



2 



solver was proposed in Refs. [T4"1[T5] . Though linearized methods are numeri- 
cally efficient and robust, they can generate unphysical states with negative 
energy density. It means that additional efforts should be made to handle 
strong enough rarefactions (see Ref. [T5]). 

HLLE method [T8][T2] is based on approximate handling of the Riemann prob- 
lem by introduction of an intermediate state and two waves propagating with 
some signal velocities. Extension of the method to relativistic hydrodynamics 
was made in Refs. [4"0fTTj . Quality of the scheme depends on choice of sig- 
nal velocities thus making it a bit of art. The scheme also does not take into 
account presence of a contact discontinuity in real solution of the Riemann 
problem. It leads to a bad resolution of contact discontinuities in particular 
and to high numerical diffusion in general. This problem is partially lifted in 
recently developed HLLC numerical scheme [4"3"ll3"4] where contact discontinu- 
ity is present in solution of the Riemann problem. HLLE method was applied 
for relativistic heavy ion collisions by Rischke et al. (38]. A higher order exten- 
sion of HLLE metod, PPM (Piecewise Parabolic Method) was introduced in 
Ref. [6] and used for description of ultra-relativistic A+A collisionsby [2"0f21~] . 

Another approach to approximate solution of the Riemann problem is two- 
shock approximation [2 7 35]. This method ignores presence of rarefactions and 
take solution of the Riemann problem consisting only of shock waves. Such 
approach works quite well for shock configurations but starts having problems 
for configurations with strong enough rarefactions due to violation of entropy 
condition. Additional work like making entropy fixes should be done to repair 
this. 

Finally, an approach that does not consider the Riemann problem at all, is 
also possible. This is flux splitting method (see e.g. Ref. [ID])- The method 
splits the flux into right going and left going parts in a way that satisfies 
conservation laws and conforms with the wave structure of the flow. One of 
the best realizations of this method is Marquina flux [33] . Marquina-type flux 
is the base of the code GENESIS [1]. 

Another problem in constructing a numerical code is the choice of accuracy. 
Accuracy is achieved by proper reconstruction of quantities inside a numerical 
cell or a proper reconstruction of fluxes if a flux splitting method is used. 
There are high order accuracy ENO (essentially non-oscillatory) [T9F 9.8] and 
WENO (weighted ENO) schemes [29123] - The recent achievements here are 
the RAM code @2] and the WHAM code [42] both of 5-th order accuracy. It 
is thought that high accuracy helps in simulating flows with different energy 
scales. 

We use a second order accuracy scheme by Falle & Komissarov [i~4"][T5] with 
a non-linear Riemann solver. This solver is constructed in such a way that 
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it is an approximate solver in case of non-zero tangential velocities and an 
exact solver in case of zero ones (including one-dimensional flows). This gives 
us an advantage of equally good treating of difficult situations such as ultra- 
relativistic flows, strong shocks and rarefactions. In our general scheme we 
incorporate the source term in right hand side of hydrodynamics equations 
aiming utilization of the code for description of the physical processes where 
such a source is required. It concerns, in particular, the matter evolution in 
ultra-relativistic nucleus- nucleus collisions. Our paper is organized as follows. 

In Section 2 the basic equations of relativistic hydrodynamics and their nu- 
merical treatment are presented. In Section 3 an exact Riemann solver for 
one dimensional flows is constructed. Section 4 contains description of treat- 
ment the special case of equation of state p (e, n) = 0. In Section 5 we present 
how our scheme interacts with vacuum. In Section 6 we give the approxi- 
mate method of including tangential velocities. Section 7 describes how to 
convert conservative variables into primitive ones. In Section 8 the method of 
including a source term is given. Section 9 is devoted to testing our code on 
several examples of one-dimensional flows. In Section 10 we test our code on 
three-dimensional flows. Section 11 summarizes the results of the paper. 



2 Basic equations 



Our task is to numerically simulate the system of equations of relativistic 
hydrodynamics. It can be written as 

<9T^ 



where T^ u = (e + p) u^u^ — pg M!/ is energy-momentum tensor of ideal fluid, 

e is energy density taken in the rest frame of a fluid element, p = p (e, n) is 

pressure, n is baryonic number density taken in the rest frame of a fluid ele- 
/ 1 v \ 

ment, u M = -, is 4- vector of fluid velocity and g^ is the metric 

V 1 — v 2 1 — v 2 ) 

tensor that in case of flat space-time has the form g M1/ = diag (1,-1,— 1,-1). 
We use such notations that Greek indices take values from to 3 and Latin 
ones from 1 to 3 (are the space coordinates). 

Equations (CQ) must be supplemented by conservation equation for baryonic 
number 

^T- = 0. (2) 
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The system (DO)-® is full and represents five equations for the set of five inde- 
pendent variables P = (e,n,v x ,v y ,v z ). This set is called primitive variables. 

The system ([I])-® is a system of conservation laws and can be rewritten as 

dQ + flF(Q) + 0G(Q) + <9H(Q) = 

dt dx dy dz 



where 



Q 



'e + pv 2 (e + p)v x (e + p)v v (e + p) 



1 - v 2 ' 1 - v 2 ' 1 - v 2 ' 1 - v 2 ' vT 



V 



2 



(Qi,Q„Q,,Q z ,Q„), (4) 



is the vector of different five independent variables that are called conservative 
variables and 

f (e + p)v x (g + P) V x (g+P)VgVy (g + P)Vo:V^ \ ( , 

{ 1-v 2 ' 1-v 2 +P ' 1-v 2 ' 1-v 2 'vT^J'^ 

/ (e + p) v y (e + p) v y Vx (g + P) (e + p) v y v z uv y \ 
u - I ~~ I ~ ' 1 — ' — Z7T~ + P' 1 3 ' / -, 2 I ' W 



1 - v 2 ' 1-v 2 '1-v 2 1-v 2 'VI 



V 



it / (g + P) v z (g + P) v ^ v x (e + p) v ^ v y (g + p) v ^ , nv ^ \ (7] 
{ 1-v 2 ' 1-v 2 ' 1-v 2 ' 1-v 2 P, VT^) ,{) 



are the vectors of fluxes in x, y and z directions respectively. 

There are numerous numerical methods for simulation of conservation sys- 
tems. As a basis we implement a Godunov-type method, described by Falle & 
Komissarov [X1][T5] in one- and two-dimensional cases. 

In order to construct a numerical scheme for system ([3]) we divide numerical 
domain into cells with size h x along the x-axis, h^ along the y-axis and h z 
along the z-axis. Knowing solution at time t n we want to find solution at time 
tn+i — tn + At. 

Integrating in time from t to t n + At and in space over domain (i — l)!^ < 
x — ih x , (j — 1)1^ < y < jh y , (k — l)h z < z < kh z we obtain the following 
exact relation 



At 



Qn+l pjn _ AAL f ™ _ rpn 
ijk ~ ^ijk i \ r ijk r (i-l)j'fc 



At / /-in 



At 
h7 



TTf 



TTt! 
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Here 



ih x jK khz 

QSfc=hhX Iff Q(x,y,z,t„)dxdydz, (9) 



Z (i-i)h a: y-i)h„ (fc-i)h a 



p 3* = Afth:/ / / F (Q(^y> z > t )) dtd y dz > ( 10 ) 

y t„ (i-i)hj (fc-i)h» 

tn+l ihj; feh z 

GS " fc = Athx/ / / G (Q( x .A»z,t)) dtdxdz, (ii) 

t„ (i-l)ha (fc-l)ll z 

t„+i ih x jh. y 

H ^ = Athx/ / / HCQ^y'*^*)) dtdxdy, (12) 



Note that in the right side of (|T0|) - (fl2|) there is integrating in time and con- 
servative numerical methods differ from each other just in ways of estimating 
these integrals. 

As a first order scheme we use the following method. We put that inside of a cell 
there is a uniform distribution of conservative quantities equal to Q^- fc . In this 
case on every cell boundary we have a discontinuity that sets up a Riemann 
problem. We have to solve it in order to estimate fluxes (TTUT) - (Tl~2T) . For this we 
solve a one- dimensional Riemann problem in every space direction. Indeed, to 
obtain, for instance, F"- fc we suppose that in a neighborhood of discontinuity 
the flow depends only on x and not on y, z, i.e. 

F «* = At / F (Q(^>AA,t))dt, (13) 



which is true as the distribution of quantities inside a cell is uniform. After we 
have solved the Riemann problem, on the place of a discontinuity we obtain 
a state Q* fQ^, Q^-i^fe) and the corresponding flux is 

F ijk = F (Q* (Qyjfc* Q(i-i)jfc)) • (14) 
Similarly we find G™- fc and H™^. 

In order to construct a second-order scheme we have to introduce a structure 
inside a cell and to let quantities change during a time interval At. 
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The second order in time is achieved by using a first order scheme to get an 
intermediate solution Q^ 1//2 at time t n + At/2 and by applying it to finding 
fluxes. 



In order to achieve the second order in space we approximate the quanti- 
ties inside a cell by linear functions. In doing so we must be cautious not 
to break scheme's monotonicity. This is important because non-monotonic 
schemes lead to artificial oscillations in neighborhood of shock waves and con- 
tact discontinuities. In order not to allow such effects an averaging function 
is used to calculate cell gradients that compares gradients in neighboring cells 
and chooses the smallest one. Doing so we have 





n+l/2 


dx j 


ijk 




n+l/2 


dyj 


ijk 


8Q\ 


n+l/2 


dz j 


ijk 



av 



av 



av 



n n+l/2 
^ijk 


n n+l/2 
^{i-l)jk 


n n+l/2 

**(i+i)jk 


n n+l/2 
^ijk 




K 


K 




n n+l/2 


n n+l/2 
^i(j-l)k 


n n+l/2 


n n+l/2 




hy 


\l y 




n n+l/2 
^ijk 


n n+l/2 
^ij(fc-l) 


n n+l/2 


n n+l/2 
^ijk 




h 2 


K 





(15) 



(16) 



(17) 



where the averaging function av(a,b) must have the following properties 



av(a,b) 



► | (a + b) , if a — > b 
0, if ab < 



a. 

b. 



if | a | / |b| — ► oo 
if |b| / |a| — > oo 



;is) 



There exist an infinite number of such functions. We use one of the simplest 
variants that has the following form 



av(a,b) 



a2 k + ?2 2 ,if a 2 + b 2 >0 
a 2 + b 2 

= 0, if ab < 



(19) 



After we have calculated gradients (|T5 l) -( fT7l) we can find second order left 
and right states for a cell boundary Riemann problem. For example, when 
estimating F™. fc we assume 



O (0 - O 

^ijk — H 



n+l/2 
ijk 



+ — 
2 



'9Q' 
dx 



n+l/2 

I 

ijk 



(20) 
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n+1/2 h x / <9Q 



n+1/2 



(21) 



(i+l)jk 



Correspondingly 



F 



<n 
ijk 



F (Q* [Q$ k , Q(ili)ifc)) • 



(22) 



Similar relations are applied to find G"- fc and H™- fc . As we said above, appli- 
cation of quantities with time index n + 1/2 in order to find fluxes warrants 
the second order of the scheme in time. 

Note that in relativistic hydrodynamics the condition 



must be held in order to obey causality principle. In general, however, applying 
relations (!20]) - (l2T]) may violate it. In such cases we have to additionally adjust 
slopes in order to hold ( 1231) . 

Thus we have a complete numerical second order scheme. The only missing 
constituent is the procedure of calculating solution of a Riemann problem i.e. 
a Riemann solver. Having it at hand we can find Q* (Q,fj k , Q^l^-^J and then 
find fluxes and then move one step further in time. 



3 Exact one-dimensional Riemann solver 

In this section we formulate a procedure of finding solution of a Riemann prob- 
lem. We call our solver exact in a sense that numerical methods are applied 
to exact relations that describe a Riemann problem solution. This differs our 
method from numerous approximate Riemann solvers that use approximate 
relations for estimating intermediate state Q*. 

As was said above, we need to construct an algorithm for a one-dimensional 
Riemann problem. For this we consider the following system 




(23) 



dQ | 0F(Q) 

dt dx 



0. 



(24) 



where 




(25) 
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wave 2 



wave 3 
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Fig. 1 . Schematic depiction of break up of an arbitrary discontinuity. Initial left Q"' and right Q' r ' states 
are separated by two intermediate states Q' and Q". All the states are connected through the set of 3 waves. 
For linear systems these waves are discontinuities. For a non-linear system of relativistic hydrodynamics 
with a thermodynamically normal equation of state wave 2 is a contact discontinuity and waves 1 and 3 
can be either shock waves or simple rarefaction waves 



(e + p) v p + ev 2 



nv 



1 - v 2 ' 1 



(26) 



As we see, all the quantities here depend only on t and x and velocity has 
only the x-th component denoted as v. The influence on solution of non-zero 
tangential velocities v y and v z will be considered later. 

For the system (12^1) we state the following initial conditions at initial time 
t = t 

, Q (0 , if x < 

Q(0,x)= ' . (27) 

Q (r) , if x > 



where ' and are some constant states. 

Quite often numerical solution of ( |2^1) - (l27j) is based on a linear approximation. 
We too base our idea on solution of a similar problem for a linear system. It is 
well known that for linear systems, i.e. systems with F (Q) = AQ, A being a 
constant matrix, solution of problem fr2^1) - fT27j) consists of 3 discontinuities (or 
n discontinuities for systems of n equations) divided by domains of constant 
flows (see Figured]). But in case of a nonlinear problem there arise continuous 
self-similar solutions called simple waves. 

At this point in order to be able to move further we need to say that only 
so-called thermodynamically normal media are the subject of present study. 
Medium is called thermodynamically normal when the following condition 
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holds 



d 2 p 



dr 2 



> 0, 



(28) 



where r = is so-called generalized specific volume and s is entropy. Con- 

n 

dition ( 1281) means that only compression shock waves and rarefaction simple 
waves are physical. Rarefaction shocks are nonphysical in the sense that en- 
tropy of a fluid element decreases when crossing a shock front. Compression 
simple waves are nonphysical in the sense that wave profile becomes singular 
at a definite moment of time. 

Media with violated condition (|28|) are called thermodynamically anoma- 
lous. The consideration here is similar to the thermodynamically normal case 
with only change that physical are rarefaction shocks and compression sim- 
ple waves. Cases when for some parameters relation (128]) holds and for some 
parameters is violated are more complicated and are the subject of further 
investigation. This is of importance because equations of state where exists a 
phase transition represent both normal and anomalous behavior at different 
values of thermodynamical parameters. 

Nevertheless from the point of view of pure mathematics both compression and 
rarefaction shocks exist as generalized solutions of equations (^^ . It means 
that we can search for pure shock wave solution of problem fr24|) -( f27|) . and at 
this point not take into account how physical are its constituents. 

Indeed, a discontinuity Q x — > Q 2 is a generalized solution of system ff2^|) if 

F(Q 1 )-F(Q 2 ) = S(Q 1 -Q 2 ), (29) 



where S is propagation speed of discontinuity and in linear case F (Q) = AQ 
is just an eigenvalue of matrix A. 

Arbitrary states and in general does not satisfy (129]) . But they can 
be sewn together through two intermediate states Q' and Q" (see Figured]). 
Then we have the following system of algebraical equations 

F (Q«) - F (QO = Si (Q« - Q') (30) 
F (Q') - F (Q") = S 2 (Q' - Q") (31) 
F(Q")-F(QW) =S 3 (Q"-Q w ) (32) 



This is system of 9 equations for 9 unknown variables - Q', Q", Si, S 2 , S 3 . 
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Once have solved it, we have a pure shock wave solution of a Riemann problem 
(1241) - (I2T|) . Taking into account that wave 2 must be a contact discontinuity 
for which p' = p", v' = v" = S 2 , the number of equations in (l3"Uj) - (13"2"j) can be 
reduced to 7. The system (1301) - (1321) is written in terms of conservative variables 
Q but in our code is solved in terms of primitive variables P = (e, v, n). In 

ultra-relativistic case we use instead of v 7-factor 7 = ^ = . So, after we 



have solved pi -(1321). we have e', n', e", n", Si, S 3 , V. 



The numerical method we use to solve (I30l) - (l3"2"l) is the Newton method. An 
obvious difficulty here is how to choose an initial approximation. As we have 
a set of 7 equations its quite a non-trivial problem. In some cases, when our 
initial approximation appears to be bad, the method fails. In such cases we 
use a switch to the alternative method of finding e', n', e", n", Si, S3, v'. 

It is based on relations that connect variables on both sides of a shock wave 
and on a well known fact that across a contact discontinuity pressure and 
velocity are continuous. First, note, that shock speeds Si, S 3 and velocity v' 
can be explicitly expressed in terms of thermodynamical variables. So, as a 
matter of fact we need 4 equations to find e', n', e", n". 

The first two equations we take from Hugoniot relations and obtain 



n' = n« 



( £ ' + p')( £ ' + p(0) 



\ (e(0+p(Q)( e (Q+p') 



n" = nM 



(£" + p')( £ " + p M) 



\ (eW+pW)( e W + p') 



Here we used the fact that on contact discontinuity p' = p". This gives us the 
third equation 

p (e', n') = p (e", n") (35) 



The fourth equation we obtain from relation v' = v". The velocity of fluid 
behind the shock front (marked by index 2) in the rest frame of fluid ahead 
of the shock front (marked by index 1) is 



y ah ( £ l) n l)£2,n2) = ± 



where the sign is defined by two factors- whether the shock is compression or 
rarefaction one and what direction does it go, right or left. 



(P2 - Pi) ( g 2 -£i) 
(£1 + p 2 ) (e 2 + Pi) 



(36) 
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(a) 



(b) 



(!) 0) "l a 
e , n , v =0 



(l) (!) n 

e , n , v =0 



Fig. 2. A shock wave g^nC'V 



•CO 



e , n , v in the rest frame of the left fluid. Here 



/ 12 (ei, ni, £2, n2) (see Eg. 1361 1'). If v^ < then we have a compression shock. When v' a > then 



there is a rarefaction wave. 

Then the fourth equation is written as 



r(0 



+ vf 2 (e(0 , n« , n') V W + vf 2 (e W , n« , e", n") 



1 + vWvf 2 (e« , n« , e', n') 1 + V M v^ 2 (eto , nM , e", n") 



(37) 



Depending on situation there can be 4 combinations of signs. 
Shocks speeds Si and S 3 can be written as 



Si 



l + v( l h- h (e( l \n®,£',n') 



(3c 



and 



vW + v+(e«,nM,e>" 
1 + v«v+ (eW,nW,e",ii") : 



(39) 



where • 



v 12 l > £i,ni,£2,n 2 ; 



(P 2 - Pi) (^2 + Pi) 



is the speed of the shock 



\[ (e 2 - ei) {si + pa) 
front in the rest frame of fluid ahead of a shock (marked by 1) and the sign 
is "-" when a shock moves in the rest frame to the left (applies to shock wave 
1) and is "+" when a shock moves to the right (applies to shock wave 3). 



Thus, we have constructed a purely shock wave solution for a Riemann prob- 
lem. But the real solution of a Riemann problem (|24|) - (|27p consists not only 
from shocks. In general, for a normal medium, waves 1 and 3 can be either 
compression shocks or rarefaction simple waves. Thus, the next step is to find 
real pattern of waves. 

As an example consider a shock transition eW,n®,v' ! ' — > e', n', v' with > n' 
(see Figure [2]). In order to find out whether such a shock is compression or 
rarefaction shock we have to go to the rest frame of the left fluid. If velocity 



12 



to the right v' a > (it means that v^ < v ; in laboratory frame) then fluid 
crosses the shock front from the left to the right. As we supposed that > n' 
it means that the shock is a rarefaction one. If v' a < (v^ > v' in laboratory 
frame) than fluid crosses shock front from the right to the left and we have a 
compression shock. Similar considerations should be applied to case n^ < n' 
and to a shock transition e^ r \ n^ r \ — > e", n", v". 

Once we have found types of shocks, our considerations are as follows. If both 
shock waves, 1 and 3, are compression shocks, then solution of a Riemann 
problem is over. If one of waves appears to be a rarefaction one, we should 
replace it by a simple rarefaction wave. In such a case we have to sew together 
a simple and a shock wave. As an example consider the case when wave 1 is 
a simple rarefaction wave and wave 3 is a shock wave. 

Profile in a simple rarefaction wave is described by the following system of 
ordinary differential equations 

dv = ± c 8 (g,n)(l-v 2 ) 
dn n 

d£ = g + P(g,n) ^ 
dn n 



where c s (e, n) 



dp 



is the speed of sound. Sign in Eq. (14131) depends on 



direction of wave propagation. If wave propagates into the left fluid than there 
is if it goes to the right fluid, there is "+". In our case of wave 1 we have 
Equations fH0|) . fHTT) define functions vf ar (n) and e rar (n). Thus, sewing of 
a simple rarefaction wave 1 and a shock wave 3 means solution of the following 
system of equations. 



(e( r ),nW,£>") +vW 
1 + v( r hf 2 (£M,nM,£",n") 



Vl2 



p(e',n') 



n 



(6" + p(e / ,n'))(e" + pM) 
\ (£(0+pM)( £ M +p ( £ ',n')) : 



(42) 

(43) 
(44) 

(45) 



where v 12 is defined by Eq. (I36[) . Thus we find e',n',e",n". Corresponding ve- 
locity is v' = v" = v mr (n'). System fT42|) - (l45|) is solved by Newton method. 
Integration of system fT4T)]) - (l4T]) to find \f ar (n) and ef ar (n) is conducted by 
using a Runge-Kutta method. 
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The velocity of wave 1 is defined by relation 



r (0 - 



Si = 1 ( 46 ) 

1 — V^'C, 



Similar considerations are applied to the case when we have a simple rar- 
efaction wave moving to the right state and a shock wave moving to the left 
state. 

In case when in all-shock solution waves 1 and 3 are rarefaction shocks, we 
have to replace both by simple rarefaction waves. In this case we solve the 
following system of equations 

v- r (n') = v+ r (n") (47) 

p(e',n')=p(e>") (48) 

e' = e rar (n') (49) 

S" = Erar (n") (50) 



Now, after we have solved a Riemann problem on a cell interface, we have 
to find state Q* (see Eq(ll4"l) and fl22|) ). In case when wave 1 and 3 are shock 
waves Q* = Q (0 when Si > 0, Q* = Q (r) when S 3 < 0, Q* = Q' when Si < 
and S 2 > 0, Q* = Q" when S 3 > and S 2 < 0. 

In case when there are simple rarefaction waves situation is more complicated. 
Take again the example with a rarefaction propagating to the left side and a 
shock propagating to the right side. Then Si is defined by Eq.( H6|) . If Si > 
then Q* = Q w . If S 3 < then Q* = Q (r) . If S 3 > and S 2 < then Q* = Q". 
But if Si < and S 2 > situation depends on how strong rarefaction is. If 
v' < c s then rarefaction is subsonic and Q* = Q'. If v' > c s then rarefaction 
is supersonic and we have to integrate system (f4"0"|) - (l4"T|) starting from the 
left-side state to the sonic point Q sn = (e sn , n sn , v sn ), i.e. to the state where 
V (n sn ) = c s . Then = Q sn . 

Similar considerations are made to find Q* in case of a shock wave propagat- 
ing to the left and a rarefaction propagating to the right and in case of two 
rarefactions. 

Now, after we have defined state Q*, construction of a one- dimensional nu- 
merical scheme is complete. In three dimensions we have to take into account 
presence of non-zero tangential velocities. This problem is considered in Sec- 
tion 4. 
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4 Special equation of state p (e, n) = 



In case of a pressureless medium construction made in Section 2 does not 
work. Here we have to apply a special algorithm. Flux we want to find is 

tn+l 



F7= At J ^ («h x ,t)) dt, (51) 

t n. 



In predictor step distribution of quantities over a cell is uniform and obviously 

F(QD, ifv?>0 
F (QIVi) , if v? +1 < ' 

In corrector step we have to take into account linear distribution of quantities 
over a cell. As an example take the case with v" > 0. Case v" +1 < can be 
taken similarly. 

In this case we take approximately 

F? = ±(f(Q?) + f(q? +1 )), (53) 

where the second term is 

F (QI l+1 ) = F (Q (ih x , t n + At)) = F (Q (ih x - Ax, t„)) (54) 

Here Ax is defined by algebraical equation 

v (ih x - Ax, t„) At - Ax = (55) 

It is solved by Newton method, v (ih x — Ax, t n ) is defined using linear struc- 
ture of quantities inside a cell. This gives us the system 

nn . ( h * A \ e (ih x - Ax, t n ) , 
Q 1 , + ai ^--Axj = 1 _ v2( .^_ AX|U (56) 

n n , n ( h x > \ _ e (ih x - Ax, t ra ) v (ih x - Ax, t ra ) 

Q " + ~ ~ l-vM,h,-Ax,t n ) ' (57) 



where a\ and a x are cell slopes of corresponding conservative variables. 
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Solving Eq. (l55l) and using Eq. (l54|) we complete construction of a scheme for 
special case of equation of state p (e, n) = 0. 



5 Boundary with vacuum 



Another difficult problem is defining fluxes at boundary of medium and vac- 
uum. Let i-th cell at time t n be vacuum and (i — l)-th cell be medium. Con- 
sider it for predictor step when distribution of quantities over (i — l)-th cell 
is uniform and equal Q£-i- Having this, we want to find Q" + . Then 



Q 



n+l 



At 



At 



Fra 
i-1) 



(5S 



where 



t n +l 

F ^ _1 = At / F (Q(( i - 1 K» t )) dt 

1 11 



(59) 



is the flux between (i — l)-th and i-th cells. In order to estimate it we use the 
following approximation. Knowing that under expansion into vacuum velocity 
must increase and all other quantities decrease and supposing that At is small 
enough, we use a linear approximation for quantities on cell boundary 



(60) 
(61) 
(62) 
(63) 



v 6 (t) 


= v?_! + a(t 


t n 


e b (t) 


= eti-P(t 


tn, 


n fe (t) 


= <i-«(t 


t n 


P 6 (t) 


= P?-i-7(t 


-t„ 



Substituting Eqs.(l6"Qll-(l6"3l into Eq.(l39l we obtain 

(eJLj + ptx) v?_ x 



K 



(l " (vf^) 2 ) 2 



2 N 



1-fv 



Av + 



(As + Ap) v?_i 



(64) 
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x(i-l) 



Pf-l+gf-l (vti) 1 
1 - (vUf + 2 



ra(i-l) 



V-l v i-l 



2v n /V n +n" 



Av+ 



Ap + Ae (v«_ x ) : 



1-fv 



+ 



v^An 



(i - (vto 2 ) 3/2 J^U) 



(65) 



(66) 



where AA = A b (t„ +1 ) — A"_ x . Take into account the following 

ih x (i-i)h !C +v uoc At 
Qr +1 = ^ / Q(x,t n+1 )dx = i- J Q(x,t n+1 )dx 

(i-l)h. x (i-l)h a 

« r^Q 1 (WO , 



(67) 



where v mc is the velocity of boundary with vacuum that at time t n +i is some- 
where inside the cell. In order to close the system we approximately put 
v mc « v 6 (t n+ i). After this we substitute ([MD-dSZD into §E§ and obtain al- 
gebraical system of 3 equations for e b (t n+ i) , n 6 (t n+ i) , v b (t n +i). After solving 
it by Newton method and using relation fl67|) we obtain Q" +1 . In a corrector 
step we use Q^i instead of Q^-i- 

The case when vacuum is located on the right side from the medium is con- 
sidered similarly. 

Note, that after we have found Q™ +1 , there will be non-zero Q2h ■ But in 
reality boundary with vacuum crosses the boundary of between the z-th and 
(i + l)-th cells in several time steps. In order to avoid such artificial superlight 
expansion of matter, we should allow matter to expand into (i + l)-th cell only 
after boundary with vacuum reaches the right boundary of i-th cell. 



6 Three-dimensional scheme: the role of tangential velocities 



In three-dimensional case for one cell we have to solve three Riemann problems 
- in x-, y- and z- directions. As an example consider a Riemann problem in 
x- direction. Then in Fig. [1] on the left and right states there are tangential 
components of velocity- and v!£\ respectively. These components 
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constitute vectors v\ and v£ . As a result, the system (I3"0"l) - (l3"2l) has 13, not 
9, equations. The structure of wave pattern is the same as in one- dimensional 
case (see Fig{T]). It is known that due to Lorentz factor coupling tangential 
velocities are incorporated into solution of the Riemann problem [3"Tf3"T] . This 
makes solution of exact equations rather numerically expensive. In order to 
simplify things and reduce to one-dimensional Riemann problem we make 
an approximation that tangential velocities are continuous over shock waves 
and simple rarefactions but undergo a jump over a contact discontinuity. It 
means that \' t = vf and vf = . Having this we can again reduce in 
( I30"l) - (l3"2l) to 9 equations and taking into account that Eq. (15111 is an identity 
to 7. To make things completely similar to one-dimensional case we must 
use Eq. (130]) in the reference frame where v[ = 0. In this case we have to 

substitute v« - vg = vW/^l-(v®) 2 and v' - v' m = v'jsjl - (vf } ) 2 . Eq. 
( 1321) thus must be used in reference frame where Vj = with substitutions 

v« - vff = ^//l-(vf») 2 , v" - vl = v:/^-^) 2 After this we 
have again one-dimensional system of equations. Correspondingly, condition 
of normal velocity continuity over contact discontinuity v' x = v" takes form 

v^i-(v«") 2 = v;; lV /i-(vr') 2 . 

All the rest of the scheme can be reduced to one-dimensional case described 
in Section 2 in similar way. For example Eq. (f3"Tj) can be rewritten as 

vff +v£(e (0 ,n ( V,n') r v^ + vt2{e^,^\e",n") r~ 

1 + v#v£ (e«, n«, e', n') V V v * ) 1+ v M y ± ( £ (r) ; nW) £ ^ n ") V 

Similarly we find fluxes in y- and z-directions and complete a time step for 
fully three-dimensional hydrodynamics equations. 



7 Conservative and primitive variables 

In Eq.flS]) conservative variables Q are updated. But in the Riemann solver 
we use primitive variables P = (e, n, v). It means that at every time step we 
have to convert conservative variables into primitive ones. 

Conservative variables are defined by Eq.(J2]) and one can express thermody- 
namical parameters e and n via velocities 

n = Q n Vl - v 2 (69) 
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^ = Qi - Q x vx - %vy - Q z Vz (70) 
Consider as an example case when Q x ^ 0. Then it is obvious that 

One can check that the following relation holds 

Qs-QiVs-p (£,11)^ = (72) 



Substituting in it relations (!69l) - (!7Tl) we have just one single algebraical equa- 
tion for v x . It can be easily solved by Newton method. 



8 Source term 



In our numerical scheme we include the possibility for presence of a source 
term. It means that we solve the following system of equations 



where S (Q, r, t) is a source term. 

In order to solve such a system we use the standard fractional step method. 
It means that as a corrector step we solve equations without source term 

d<y + dF(Q p ) + dG(Qn + dH(Q p ) = q 

dt dx dy dz 



We just have constructed the procedure how to solve this system. As a result 
we obtain (Q p )™ +1 - Then we use these values as initial conditions for the 
following system of ordinary differential equations 

^ = S(Q,r,t) (75) 



To solve the system (1751) we use a Runge-Kutta method. Time step is At. This 
solution gives us numerical solution of system (1751) . 
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9 Test simulations: One dimension 



9.1 Shock tube problems 



Here we consider two shock tube problems. Equation of state is p = c 2 c e with 
cl = 1/3. 

The first problem is collision of two gas masses that gives two shock waves 
propagating to opposite sides. To model this situation we take the following 
parameters: e {l) = 2, vW = 0.7, = 2 to the left and e (r) = 4, v^ = -0.5, 
n( r ) = 4 to the right. Results of numerical simulations are depicted on Fig [31 

The second problem is the contact of two motionless masses of gas with dif- 
ferent densities. This can give a rarefaction wave and a shock wave. To model 
this situation we take the following parameters: = 10, v^ = 0, = 10 
to the left and e^) = 2, = 0, n^) = 2 to the right. Results of numerical 
simulations are depicted on FigHl 



9.2 Bjorken-like flow 



In this subsection we numerically simulate Bjorken-like expansion solution. 
This solution was introduced and analyzed in Refs. [22"|T7|3] . Equation of 
state is p = c 2 s e. Solution has the following form 



so 



T 



T 



(76) 



n 



n - 



r 



T 



X 

V= ? 



(77) 
(78) 



where r = 




We take the following parameters: c^ = 1/3, r = 5, n = 1 and e = 1. Initial 
time is to = 5 and final time is t = 8. At moment of time to we take all 
the quantities in region |x| < t — a, where a is small enough, being defined 
by relations (1761) - (1751) . Outside this region all the quantities are put to zero. 
Results of numerical simulations are depicted on Fig. [5j 
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(a) 



— analytical curve 
numerical curve 



-analytical curve 
numerical curve 



(b) 



(d) 



- analytical curve 
numerical curve 



2 4 6 



analytical curve 

■ numerical curve 



2 4 



-analytical curve 
numerical curve 



(f) 



analytical curve 

■ numerical curve 



Fig. 3. Results of numerical simulations for a shock tube problem with following parameters: eW = 2, 
vW = 0.7, nW = 2, £< r ' = 4, v< r ) = -0.5, n< r ' = 4. Initial moment of time is to = 0, final point of time 
is t = 3. All quantities are taken at time t. (a)-(c) are energy density, velocity and baryonic charge density 
correspondingly, (d)-(f) arc the same quantities but calculated with better spatial resolution, (a)-(c) have 
100 grid points and (d)-(f) have 1000 grid points in space. 
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numerical curve 



- 



(b) 
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0,30- 
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(c) 



- analytical curve 
numerical curve 




(d) 
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numerical curve 
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(e) 

0,35 
0,30 
0,25 
0,20 
V 0,15 
0,10 
0,05- 
0,00 
-0,05 
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numerical curve 



(f) 



— analytical curve 
numerical curve 



4 6 8 10 



Fig. 4. Results of numerical simulations for a shock tube problem with following parameters: gW = 10, 
v (i) = 0i n (i) = 10i £ (r) = 2j v (r) = o, n< r ) = 2. Initial moment of time is to = 0, final point of time is 
t = 5. All quantities are taken at time t. (a)-(c) are energy density, velocity and baryonic charge density 
correspondingly, (d)-(f) are the same quantities but calculated with better spatial resolution, (a)-(c) have 
100 grid points and (d)-(f) have 1000 grid points in space. 
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(a) (b) 

analytical curve analytical curve 




analytical curve 

■ numerical curve 




o.o- -Vm 

-10 -5 5 10 

X 

Fig. 5. Results of numerical simulations for a Bjorken-like flow, (a) is energy density, (b) is velocity and (c) 
is baryonic charge density. Evolution starts at time to = 5 and is finished at time t = 8. All the graphs are 
taken at time t. One can see a rarefaction wave expanding into vacuum and distorting analytical solution 
on the edges. In the middle of the region numerical and analytical solutions coincide. 

9.3 Hubble-like flow 



Another example is Hubble-like flow. It is constructed similarly to Bjorken-like 
solution, see Ref. [5], except that Hubble-like solution is spherically symmet- 
ric. It means that x in system fl24|) means radial coordinate and here arises 
geometrical source term S (Q,x, t) that has the following form 

snx , t) = UJ^,-^^,-^L_). (79) 

V x(l-v 2 ) x(l-v 2 ) xVl-v 2 / 



Solution itself has the following form 

e = e [^y (80) 

n = n (- 
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X 



Fig. 6. Results of numerical simulations for a Hubble-like flow, (a) is energy density, (b) is velocity and (c) 
is baryonic charge density. Evolution starts at time to = 5 and is finished at time t = 8. All the graphs are 
taken at time t. One can see a rarefaction wave expanding into vacuum and distorting analytical solution 
on the edges. In the middle of the region numerical and analytical solutions coincide. The difference is at 
x = because source term l|79|l has there singularity. 



(82) 



where r = y t 2 — x 2 . 

We take the following parameters: c 2 = 1/3, r = 5, n = 1 and e = 1. Initial 
time is t = 5 and final time is t = 8. At moment of time to we take all the 
quantities in region |x| < t — a, where a is small enough, being defined by 
relations (!80l) - (l82l) . Outside this region all the quantities are put to zero. An 
artificial extension of solution is made into region x < in order to avoid 
setting boundary conditions at x = 0. Results of numerical simulations are 
depicted on Fig. El 
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9.4 Spherically symmetric solution with p = 



In the paper by Sinyukov & Karpenko [H] a new class of 3D anisotropic 
analytical solutions with quasi-inertial flows has been found. Here we consider 
a spherically symmetric solution for expanding finite system with equation of 
state p = that belongs to this new class. For testing of our code we choose 
the following specific form of the solution 



\(t+t x y 



(t + t, 



(t+t.) 



(83) 



x 



t + t a 



J4) 



n 



(t + t a 



b^t+t,) 
'(t+t«) 



-X 2 



Here as in previous subsection exists geometrical source term. We take the 
following parameters: c e = 1500, c n = 1500, h 2 e = 0.5, h 2 n = 0.5 and t x . = I. 
Initial time is to = and final time is t = 3. At moment of time t we take 
all the quantities in region |x| < t + t x — a, where a is small enough, being 
defined by relations (jS3]) - (j81|) . Outside this region all the quantities are put to 
zero. An artificial extension of solution is made into region x < in order to 
avoid setting boundary conditions at x = 0. Results of numerical simulations 
are depicted on Fig. [7J 



10 Test simulations: three dimensions 



10.1 Hubble-like flow 



Here we again use as a Hubble-like flow (|80p - (|82p as a test. But, unlike the 
previous section, we use fully three-dimensional code to simulate it. Of course, 
for three-dimensional calculations there is no source term. The results of nu- 
merical calculations are depicted on Figs. [HJ19J Parameters of flow are taken the 
same as in subsection 9.2. We show dependences on distance from center along 
x-axis, in xy-plane for x = y and in a radial direction for which x = y = z. 
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X X 

Fig. 7. Results of numerical simulations for a spherically symmetric flow with p = 0. (a) is energy density, 
(b) is velocity and (c) is baryonic charge density. Evolution starts at time to = and is finished at time 
t = 3. All the graphs are taken at time t. In graphs (a)-(c) number of cells in space grid is 100. In graph (d) is 
depicted energy density simulated with number of cells in space grid being 1000. It is seen that irregularities 
that exist for 100 grid cells almost disappear if take 1000 grid points 

10.2 Elliptical solution for p = 



Here we take a solution similar to that taken in subsection 9.4. The difference 
is that now in order to test three-dimensional code solution is not spherically 
symmetric. In accordance with results of [JT] we choose the solution in the the 
following form 

P. be 

e~^s (86) 



(t + t x ) (t + y (t + 1,) 



X 



t + tj t + tj,' t + t a 



n = Cn 

(t + t,) (t + t„) (t + t z ) 
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Fig. 8. Results of numerical simulations for a Hubble-like flow in three dimensions, (a)-(c) are energy 
density, velocity and baryonic charge density respectively along x-axis. (d)-(f) are the same quantities in 
xy-plane for x = y. (g)-(h) are energy density and velocity in a radial direction for which x = y = z. 
Baryonic charge density is depicted on Fig [9] Evolution starts at time to = 5 and is finished at time t = 8. 
All the graphs are taken at time t. One can see a rarefaction wave expanding into vacuum and distorting 
analytical solution on the edges. In the middle of the region numerical and analytical solutions coincide. 
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analytical curve 
numerical curve 



Fig. 9. Baryonic charge density in a radial direction for which x : 
mensional simulations. 



z for a Hubble-like flow. Three-di- 



where s 



X 



t + t, 



t + t. 



t + t 



This is a generalization of 



spherically symmetric flow considered in subsection 9.2. 



The results of numerical simulations are depicted on Figs [TOinTl We take the 

1.5 and 



1500, 



1500, V e 



0.5, hi 



0.5, ty 



following parameters: c e = 
t z = 2. Initial time is to = and final time is t = 3. At moment of time to we 
take all the quantities in region |x| < to + t z — a, where a is small enough, 
being defined by relations fl86l) - fl88l) . Outside this region all the quantities are 
put to zero. 



10.3 Hubble-like flow with a source term 



Here we consider solution of the following form 



e = e 



'to (r,t) 



3(1+C2) 



39) 



n = n 



r 



V (r,t)' 



(90) 
(91) 



where r = \/t 2 — r 2 . Equation of state is p = cje 



28 




Fig. 10. Results of numerical simulations for an elliptical pressureless flow, (a)-(c) are energy density, ve- 
locity and baryonic charge density along x-axis respectively, (d)-(f) are the same quantities along y-direction 
and (g)-(h) are energy density and velocity along z-direction. Baryonic charge is shown on Fig llll Evolution 
starts at time to = and is finished at time t = 3. All the graphs are shown at time t. 
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Fig. 11. Results of numerical simulations for an elliptical pressurelcss flow, (i) is baryonic charge density 
along z-dircction. (j)-(l) are energy density, velocity and baryonic charge density in xy-plane for the line 
x = yrespectively. (m)-(o) are the same quantities for a radial direction for which x = y = z. Evolution 
starts at time to = and is finished at time t = 3. All the graphs are shown at time t. 
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The fact that To is not a constant leads us to an artificial source term S 
(S e , S x , S y , S z , S n ) where 



3(1 



c 2 a )et 2 



T T Z 



1 + 



.2 



dr 
dt 



l + c; 



t dxj 



(92) 



S; 



3(1 + cg)et 
t t 



3nt 



dr 
dt 



l + cf)x» dr (l + c 2 s )x l • dr r «9r, 



T 



X 1 8t 

7 



at 



+ 



T t 



- + 



Stf ' t dxi 
i = x,y,z 



(93) 
(94) 



We take the following parameters: c^ = 1/3, no = 1 and £o — 1- Initial time is 
to = 5 and final time is t = 8. At moment of time to we take all the quantities 
in region |r] < to — a, where a is small enough, being defined by relations 
(I9"2l - (I9"1"|) . Outside this region all the quantities and source term are put to 
zero. Function r (r, t) is taken to be (a) r (r, t) = t and (b) r (r, t) = e~ axy 
with a = 0.1. Results of numerical simulations are depicted for case (a) on 
Fig. [12] and for case (b) on Fig. 



10. 4 Source term defined by an arbitrary flow 



Suppose that we have an arbitrary hydrodynamical flow. It is characterized 
by profiles of energy density e(r,t), velocity v (r, t), baryonic charge density 
h (r, t) and by equation of state p (e, n). Then this flow is a solution of hydro- 
dynamics equations with a source term 



dT 



dx u 



S„ 



(95) 
(96) 



where 
S" 



dx u 



<9nu M 



(97) 



and = (e + p^u" - p 
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Fig. 12. Results of numerical simulations for a Hubble-like flow with a source term defined by to = t (see 
Eq J89l l- H91l l). (b)-(d) are energy density, velocity and baryonic charge density along x-direction. (e)-(g) are 
the same quantities in xy-plane for the line x = y. Evolution starts at time to = 5 and is finished at time 
t = 8. All the graphs are taken at time t. (a) are initial values of energy density at time to = 5. Outside the 
region with matter there is no source. It can be seen on (b)-(d) that inside the region with matter numerical 
and analytical solutions coincide but outside this region solution is distorted. Namely, expansion there is 
faster that in inner region as choice of to = t makes expansion in inner regions slower. 
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Fig. 13. Results of numerical simulations for a Hubble-like flow with a source term defined by to = e~ axy 
(see Eq. (S9J - (9TJ ) . (a)-(c) are energy density, velocity and baryonic charge density along x-direction. (d)-(f) 
are the same quantities in xy-plane for the line x = y. Evolution starts at time to = 5 and is finished at time 
t = 8. All the graphs are taken at time t. Note that solution is extremely asymmetric but numerical solution 
gives good agreement with analytical one in the region where source term is non-zero (see explanations for 
Fig. [12}. 
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One of the simplest choices of a flow is 



e = e e ~r~ 



(98) 



r 

v 



(99) 



t 




(100) 



where r = Vt 2 — r 2 . For such a flow one can obtain 



^4 + ^ 

T \ T 



) 



(101) 



s 



n 




) 



(102) 



For numerical simulation we took the following parameters- equation of state 
p = cj.e, c 2 = 1/3, b e = 1, e = 100, b n = 1, n = 100. Initial time is t = 5 
and final time is t = 8. At moment of time to we take all the quantities in 
region |r| < t — a, where a is small enough, being defined by relations (198]) - 
(jlOOj) . Outside this region all the quantities and source term are put to zero. 
Results of numerical simulations are depicted on Fig [TH 



11 Conclusions 

We presented a numerical code with a new Riemann solver. This solver is exact 
for one- dimensional flows. Exactness means that we use real wave patterns in 
order to obtain fluxes. Such an approach allows us to treat correctly vast range 
of velocities and configurations. Testing shows that the scheme works correctly 
in ultrarelativistic regimes and is capable of correct catching strong shock 
waves as well as strong rarefactions. In this sense our scheme is universal and 
can be applied to simulations of processes with large gradients of quantities 
and with large difference in velocities. Also, it account for possible source 
terms in right hand side of hydrodynamics equations. And the fact that the 
method uses physical considerations makes it especially reliable. On the other 
hand due to necessity of solving non-linear algebraical equations our code is 
more numerically expensive than the methods like HLLE or flux splitting even 
though no matrix decomposition is needed. But drastically increased efficiency 
of computers makes it possible to apply methods like ours. 
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(a) (b) 




Fig. 14. Results of numerical simulations for an arbitrary flow defined source term (see Eq. l|98| l- lll00| l). 
(a)-(c) are energy density, velocity and baryonic charge density along x-direction. (d)-(f) are the same 
quantities in a radial direction for which x = y = z. Evolution starts at time to = 5 and is finished at time 
t = 8. All the graphs are shown at time t. On (a) we depict also initial values of energy density (lower solid 
line graph). It is done in order to show the region of non-zero source term (it is non-zero in region where 
matter is non-zero at initial time). It is seen that in this region analytical and numerical solutions coincide. 
Outside this region analytical source term solution is distorted as there is no source term. 
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